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Abstract 

The radial return and Mendelson methods for integrating the equations of classical plasticity, which 
appear independently in the literature, are shown to be identical. Both methods are presented in detail as 
are the specifics of their algorithmic implementation. Results illustrate the methods’ equivalence across a 
range of conditions and address the question of when the methods require iteration in order for the plastic 
state to remain on the yield surface. FORTRAN code implementations of the radial return and Mendelson 
methods are provided in the appendix. 


1. Introduction 

With the advent of modem computers and commercial finite element analysis (FEA) codes, the 
inelastic analysis of structures is economical and accessible to design engineers. Accounting for 
inelasticity is critical in today’s simulation-based design paradigm in order to capture the true nature of 
the structural response. This is particularly true in the case of aerospace structures where weight is the 
driving design parameter. Allowing local inelastic deformation of materials provided the effects of this 
deformation can be accurately modeled and shown to be safe, can increase design efficiency by 
minimizing overly conservative designs intended to stave off yielding. As such, mathematical models that 
treat inelastic deformation of materials have become increasingly important over the last 50 years. 
Obviously, metallic materials are the primary application of such models. However, polymers and other 
non-metallic materials that exhibit permanent deformation are beneficiaries of plasticity models as well. 

The basic structural problem of interest is the boundary value problem consisting of a specified 
structural geometry along with specified loading acting on the structure. The goal is to determine the 
stress, deformation, and strain distribution throughout the structure, which is accomplished through 
solution of the equilibrium equations while enforcing the boundary conditions, compatibility, and the 
constitutive relations for the material. In the general, nonlinear case, initial conditions must be specified 
and, because the problem may be time/history dependent, the solution must be accomplished through 
incremental and/or iterative methods. As mentioned above, FEA codes provide a convenient means to 
solving such problems, although many useful structural problems can still be addressed via conventional 
analytical methods. 

This paper addresses one part of the above posed general boundary value problem; namely, the 
constitutive relations for the material. Equally applicable to FEA and conventional analytical solutions, it 
will be assumed herein that the constitutive relations in question are active at a point within a structure 
and that the structure-scale solution satisfies the equilibrium equations, compatibility, and the specified 
boundary conditions. The problem of interest thus reduces to: 
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• Given an admissible specified stress/strain state at a given location 

• Determine the unspecified stress/strain components at that location 

For instance, in the case of FEA, it is the role of the finite element method to ensure that the equilibrium 
and compatibility equations are satisfied. At a given element integration point and a given loading step, n, 
for a given strain increment, As, the problem of interest becomes determination of the new stress state at 
the integration point for loading step n+ 1. This obviously assumes that the converged state is known 
completely at step n. 

We now further restrict ourselves to the quasi-static (i.e., non-dynamic) and small deformation 
regimes and concentrate on rate-independent classical plasticity. Classical plasticity is popular due to its 
simplicity and ease of characterization. In the case of isotropic hardening, a single stress-strain curve is all 
that is required to characterize the classical plasticity parameters. Obviously, this simplicity limits its 
accuracy, but classical plasticity does capture the first-order effects of the plastic phenomenon and is quite 
useful in many instances. Classical plasticity employs a yield function, which describes the onset of 
plasticity, and an associated flow rule, which dictates the evolution of the plastic strain components. In 
some cases, local iterations are necessary in order to ensure that the stress state remains on the yield 
surface during plastic deformation. These iterations are independent of any global iterations required to 
solve the structural problem and are needed so that the correct local stress/strain state (at a point in the 
structure) are determined. 

The algorithmic implementation of the classical plasticity equations, their integration, and the need 
for the above-mentioned local iterations are the main focus of this paper. We consider two 
sources/methods that describe classical plasticity equation integration. The first is the monograph by 
Mendelson (1968, 1983). Mendelson’s integration method has been employed by the present authors for a 
number of years in a series of non-FEA based micromechanics theories and codes for the inelastic 
response of composite materials. The second is the so-called “radial return method” (Simo and Taylor, 
1985; Simo and Hughes, 1998), which appears to be much more popular and well-known than 
Mendelson’s method. The motivation for this investigation was to determine whether some advantage 
(computational or otherwise) could be gained by employing the radial return method rather than 
Mendelson’s method, given the popularity of the former. As will be shown, such an advantage was not 
found because the methods are, in fact, equivalent. While it is not the intent of this paper to attempt to 
discern which method arose first, we note that the roots of the radial return method are generally 
attributed to Wilkins (1964). However, it appears that both methods were developed independently. 

Restricting the discussion to isotropic hardening (both linear and nonlinear), we begin by presenting 
in detail the radial return and Mendelson methods, including details of the integration algorithms. Next 
the two methods are explicitly shown to be equivalent. Note that this equivalence is summarized in 
tables 1 and 2 via side by side comparisons of the steps involved in the two methods. Results applicable 
to a monolithic (unreinforced) elastoplastic material are then presented for the two methods to indicate the 
methods’ equivalence while investigating the role of the local iteration procedure. Finally, results are 
presented for a composite material whose matrix material is treated as elastoplastic within a 
micromechanics model developed by the present authors. It should be noted that explicit nonlinear 
isotropic hardening, which was included in the radial return method by Simo and Taylor (1985), has now 
been included within the Mendelson method. Mendelson’s (1968, 1983) hardening is restricted to the 
linear and piece-wise linear cases. In addition, the FORTRAN code for the radial return and Mendelson 
methods employed for the presented monolithic results is given in the appendix. This code, with some 
modification, is suitable for incorporation within an elastoplastic structural analysis code. 
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TABLE 1.— STEP BY STEP SUMMARY COMPARISON OF THE RADIAL RETURN AND MENDELSON METHODS 


Step 

Radial return method 

Mendelson method 

1 

Compute elastic trial stress from equation (6), 
s »+i = s „ +2G(e„ +1 -e„) 

Compute elastic trial strain deviator (referred to as the 
modified total strain deviator by Mendelson), e' , from 
equation (32), and the equivalent modified total strain, £et , 
from equation (37), 

e=e n+l _E n P £ et = Jj e !/ e !j 

2 

Check the yield criterion, equation (10). 

/( s« +1 .K) = S^ +1 - ^k(e P 'J = 0 

If /( sJ +1 ,K) < 0 , set s„ +1 = S J +1 , £ p +l = £ P , and the increment 
is completed. Otherwise, continue. 

Check the yield criterion, which can be written as, 

/(A K ) = A< ^ ) = 0 

If /(e',v) < 0 , set s„ +1 =2 Ge’, £ p +l = £ p , and the increment 
is completed. Otherwise, continue. 

3 

Calculate the unit normal vector, 

I /II T II 

n = s »+i/|Vn|| 

and determine the value of yAf from equation (17), 

s »+i \\-2GyAt- k(?^ +[ ) = 0 

In the case of linear hardening, the result is, 

\\< + l -\/2/3( 7 + fff ,J EE t 

Calculate the modified proportionality constant, AT' , from 
equation (52), 

au = i- k (^') 

3Gs e , 

In the case of linear hardening, the result is, 

1 I" Y + HE 1 

AU = 1 p " 

l + H/3G[ 3 Ge e , J 

see table 2 for the nonlinear hardening case. 

7 2G + 2/3H E-E t 

see table 2 for the nonlinear hardening case. 

4 

Compute the equivalent plastic strain at increment n + 1 
from equation (6) as follows, 

E Pml =E P ,„ +V2/3yA/ ( 18 ) 

and the plastic strain components from equation (5) as 
follows, 

<+i =e£ +yAtn (19) 

Compute the equivalent plastic strain at increment n + 1 
from equation (41), 

£ p,„ + , = £ Pn + AX'E el 

and the plastic strain components from equation (43) as 
follows, 

< + i= E » +AUe' 

5 

If the hardening is nonlinear, check for global 
convergence by comparing the values of | Sjj+] | detennined 

from the yield function (1) and from the constitutive 
equation (7), 

H-iH 5 -.,.) < 20) 

If the hardening is nonlinear, check for global 
convergence by comparing the value of g n+1 detennined 

from the yield function, o„ +1 i ) = 0 , to the value of 

a (I+1 detennined from the constitutive 
equation (7), 

1 ~ 2G^j—^ej /1+ j ~ £jj f n+l'j( e i/,n+l ~ 

If (jff+i - aj ^.1 > TOL go to step 2. Note that in the case of 
linear hardening with loading the is not completely 
applied strains, an additional convergence criterion is 
needed to ensure that the applied s„ +1 is equal to a Y n F +i and 

-c 

CT »+ i • 

| s n+l | “ 2G^ejj' n+ i — ^ ( y_„+i — & ij,n+ lj 

If | s -| - s c +1 > tol go to step 3. Note that in the case of 

linear hardening with loading that is not completely 
applied strains (i.e., blended loading), an additional 
convergence criterion is needed to ensure that the applied 
s„+i is equal to s ^, and |£ +1 ||. 

6 

If necessary, calculate the deviatoric stress from 
equation ( 1 1 ) or equation (4), 

s „+i = Jj K (u>,„ +1 )" or s »+i =2G ( e *+i- £ «+i) 

If necessary, calculate the deviatoric stresses from 
equation (68) or equation (7), 

or s -i= 2 G ( e » + .-< +1 ) 

J b et 
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TABLE 2.— STEP BY STEP SUMMARY COMPARISON OF THE RADIAL RETURN AND MENDELSON METHOD 
PROCEDURE FOR DETERMINING yAt AND AX' IN THE CASE OF NONLINEAR HARDENING 

NEEDED IN STEP 4 OF TABLE 1. 


Step 

Radial return method 

Mendelson method 

1 

Calculate for iteration k. 

Calculate for iteration k, 

E p*n + 1 =E p,n+ 1 + (AU) S ef 

2 

Calculate the derivative of the function g (with respect to 
yAt ), as follows (see section 2.2 for details), 

g(y A 0= |^+i -2Gy A ?- 

1 / 2/3 [k ( z Pn ) + V 2/3 yAt k' ( e p J] 

g ,= -| k '(u>J- 2G 

Calculate the derivative of the function g (with respect to 
AX' ), as follows (see section 4.2 for details), 

g(AU) = G-— 1 — [k^J + AUs^k'^JI-GAV 

3 

By the Newton-Raphson method, 

. ,*+1 , , k g[(y A 0 A ] 

(yAt) — (yAt) L J 

g'[(y A d J 

By the Newton-Raphson method, 

(M,) -<“■! -..[(V] 

4 

!f g[~(yAt)* 1 > TOL then k-*k+\, go to (1) 

If g|~(Auf 1 >TOL then k — > k+\, go to (1) 


2. Radial Return Method 

Let us define the von Mises yield criterion with isotropic strain hardening as follows, 


/( s,k) 



= 0 


where s is the deviator of the stress a, namely, 

s ij ~ ®ij ~ 


( 1 ) 


( 2 ) 


and is the Kronecker delta. Note that for convenience and clarity, a mix of vector and indicial notation 
will be employed. In equation (1), the norm of s is given by, 

||s|| = 7^ = ^ (3) 

where J 2 is the second invariant of s. In addition, j is the isotropic strain hardening rule, which 
depends on the equivalent plastic strain: 

^ = (4) 


where s jj are the plastic strain components and A indicates an increment, which may accumulate between 
step n and step n + 1 of the imposed loading profile. The flow rule is given by, 
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(5) 


As fj - yA t 


df_ 

ds ij 


- yAtn 


where n is the unit vector, n = s/||s|| , and yA t is the magnitude of the plastic strain increment, which is 
in keeping with the nomenclature of Simo and Taylor (1985). 

Substituting equation (5) into equation (4) yields. 


As, 



(6) 


The deviatoric nature of the plastic strain allows the constitutive equation for isotropic elastoplastic 
materials to be written as, 


s = 2G(e-e p ) 


(V) 


where G is the shear modulus and e is the deviator of the total strain s, 

e ij = S ij~^ E kk^ij 


( 8 ) 


According to the radial return mapping algorithm, at the current increment, n, we assume an elastic trial 
stress applicable to the next increment, n+ 1, 

sLi = s„+2G(e„ +1 -e„) (9) 


where e„ +1 is the deviatoric strain at the next increment. This deviatoric strain may be known exactly, as 
in the case of pure strain specified loading on the material, it may be derived from blended stress/strain 
loading on the material, or it may be passed to the plasticity model from a higher scale model such as the 
finite element method. The radial return algorithm requires that the correct e ;!+1 is known or can be 
determined through a global iteration procedure. 

From the yield criterion, equation (1), if 


/(«» + 1>K) 




( 10 ) 


T 

then no yielding occurs (i.e., the deformation is elastic) and s„ +1 = s „ +1 . Otherwise, from equations (5) 
and (7), we have, 


S»+l _ 2 q I *b;+l £ n _ " 

At l At 

Substituting for s„ in equation (11) using equation (9) we get, 


( 11 ) 
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s„ +1 =s^ +1 -2GyAtn 


( 12 ) 


Since n = s/ s , we have, s 


»+i 


= s 


n + 1 


n , and by substitution into equation (12), it can be shown that the 


unit normal vector n can be determined in terms of the trial elastic stress s n+ \ according to, 


J 

s n+ 1/ 


s «+l 


. This fact is key to the radial return method because it indicates that the directionality of 


the trial stress is identical to that of the converged stress, which is on the yield surface. It is therefore 
possible to return to the yield surface along this unit vector. 

By multiplying equation (1) at increment n + 1 by n , we obtain, 




n = 0 


(13) 


Thus, from s„ +1 = s„ +1 n , it follows that, 


s„ + 1 -J-k(s 


p ,„+ i 


n = 0 


(14) 


From equation (12), 


»n+ 1 


S H + 1 


-2GyAth = n 


k t 

s »+l 


-2GyAt 


(15) 


Substituting equation (15) into equation (14), we obtain, 



2GyAt- 




(16) 


Flence, we will define the function g(yAt) as the bracketed term in equation (16), namely, 


g(yAt) 



2GyAt 



(17) 


This nonlinear equation is used to determine the value of yA t . 

With the above equations, the formulation required to determine all variables at increment n + 1 given 
their values at increment n has been completed. The radial return method can thus be summarized as 
follows (see also table 1): 

(1) Compute elastic trial stress from equation (6), 

sLi = s„+2G(e„ +1 -e„) 


(2) Check the yield criterion, equation (10). If /(s^ + 1 ,k) < 0, set s„ +1 = s^ +1 , £ ; ^ +1 = , and the 

increment is completed. Otherwise, continue. 
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(3) Calculate the unit normal vector, 


T T 
n — s ;j+l/ s h+1 


and determine the value of yA t from equation (17). See sections 2.1 and 2.2 below for the linear and 
nonlinear hardening cases. 

(4) Compute the equivalent plastic strain at increment n + 1 from equation (6) as follows, 


& Pv,+l ~ 8 Pm + "J 3" YA? 


and the plastic strain components from equation (5) as follows, 

< + i=<+yAtn 

(5) If the hardening is nonlinear, check for global convergence by comparing the value of ||s„ +1 || 
determined from the yield function, equation (1), to the value of ||s„ +1 | determined from the 
constitutive equation (7). That is, the value of ||s„ +1 II from the yield function is given by, 


K ( 8 /Wi) 


and the value of s„ +1 from the constitutive equation (7) (utilizing eq. (3)) is, 


^n+\ £y 9 n + 1 {^ij,n + 1 ^ij,n + 1 


YF j. 

If s n+1 - s„ +1 > TOL go to step 3. Note that in the case of linear hardening with loading involving 

anything other than six specified strain components, an additional convergence criterion is needed to 

YF C 

ensure that the applied s, J+1 is equal to s„ +1 and s„ +1 . 

(6) If necessary, calculate the deviatoric stress from equation (14) or equation (7), 


2 l- 


S„_ L i = , - K 


( S P„, + l) A 


or s„ +1 =2 G e„ +1 -s^ +1 


2.1 Linear Hardening 


In the case of linear hardening, in which, 


' It is suggested to use a TOL value that is a small fraction of s„ +1 , such as TOL = 0.0000001 s n+ 
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k (s p ) = 7 + //8 / 


( 22 ) 


where 7 and Hare the yield stress and hardening slope (i.e., H - da/8s p , where a = ^3/2 .v // .v ;/ is the 
effective stress) in simple tension, we can readily determine the value of yAt from equation (17). Because 
equation (22) is linear, the function k at step n + 1 is known exactly in terms of the value of k at step n. 


K 


( 2 * * * 6 /Wi) = K ( e /J 


+ H 


8 _ - 8 ,, 
/^w+l Pv 


(23) 


Substituting equations (22) and (18) into equation (23) yields. 


K 


(°p,J= y+h ‘p., +h s! 


(24) 


and substituting equation (24) into equation (17) gives, 


yA t ■ 


S 72 + l 




Y + H&, 


2G + -H 
3 


(25) 


Note that the relationship between H and the actual slope (the tangent modulus, E T ) of the uniaxial stress- 
strain curve in the plastic region is given by, 

(26) 

hj — hj 1 

where E is the elastic tensile modulus. Linear hardening can also be implemented as piece-wise linear 
hardening in order to approximate nonlinear hardening. 


2.2 Nonlinear Hardening 

In the general case of nonlinear hardening the value of yA t is obtained from equation (17) by 
employing the Newton-Raphson algorithm as follows (see also table 2): 

(1) Calculate for iteration k, 


—k+\ _ —k+\ , 

& P , n + 1 “ 8 P,«+ 1 + 


Vf (yA;) ‘ 


(2) Calculate the derivative of the function g (with respect to yA t ), as follows, 

From equation (6), 


1 P 


*v... +J 


Expanding the function k as a lst-order Taylor series, 
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Hence, from equation (17), 


g(yAf) = 
Therefore, 


*n + 1 


-2GyAt- J- 


K 


( s pJ + \§~' A,K '( i 




(3) By the Newton-Raphson method. 


-2G 


{yAt) k+l = (yAtf -- 


g 

1 1 

< 

1 1 

g' 

1 1 

<1 

1 1 


(4) If 




(yAt) A > TOL 1 then &— >■ k+\, go to (1) 


Let us illustrate the nonlinear hardening for the case in which the hardening function is given by, 


K 



Y 


H_ 

A 



(27) 


where Y is the yield stress in simple tension and H and A are material parameters. Elastic-perfectly plastic 
material response is obtained when H- 0 and bilinear material response (i.e., linear hardening) is 
obtained when A is small (but not zero), yielding 

k(s p ) = Y + Hs p (28) 

and we see that the parameter H is then the hardening slope. In the general case, H is the initial slope of 
the nonlinear portion of the material response. Finally, for large values of s p , the stress will saturate to a 

value of Y + HI A. Note that the derivative, k'(s /)h j , needed in step 2 of the radial return procedure 
above, is given by, 

K'(e p J = //e-^ (29) 


* A small value for the Newton-Raphson tolerance, such as 0.0001, is suggested. 
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3. Mendelson Method 


Mendelson (1968) derived plastic strain - total strain plasticity relations, which enable the 
computation of the plastic strain increments from total strains without recourse to the stresses. As in the 
radial return method, it is assumed that the load is incremented from a completely known state at 
increment n to the next increment n+ 1, with the task being to determine the state at this next increment in 
the context of classical plasticity. The strain at the next increment can thus be written as, 

e„+i =e e n+l +e% + At p (30) 

where the superscript “e” refers to an elastic quantity. Subtracting the mean strain ( & mean = s u , /3 ) from 
both sides of equation (30) yields, 


e„ +1 =< +1 +££ + A z p (31) 

The Mendelson method defines the modified total strain deviator associated with the increment as, 


Combining equations (32) and (31), 


p f — p p P 

c ^n + 1 c n 


e' = ef i+1 +A£ p 


(32) 


(33) 


From Hooke’s law and the Prandtl-Reuss equations we have, 


®«+l 


Az p 


2 G 2GAX 


»n+l 


(34) 


where A/, is the Prandtl-Reuss proportionality constant. Substituting equation (34) into 
equation (33), 


e' = fl + — - — 4 Az p 
V 2GAXJ 


(35) 


Multiplying both sides of equation (3 1 ) by 2/3 times itself, 


2 , , 2 

— P-P - = — 

3 lJ lJ 3 


1 + - 


1 


2GAX 


Ae P jAs p . 


(36) 


Mendelson defines the equivalent modified total strain as, 


'et 


c r - 

V u 


( 37 ) 


which, along with equation (4) allows the following equation to be written from equation (37), 
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J et 


i+- 


\ 

2GAA, 


As, 


We now define a modified proportional constant AA' such that, 


1 =. + - 1 


AA' 2GAA 


so equation (39) can be written as, 


Then, 


Az p =AVs et 


s a„ + i= s a„ +A ^ 


(38) 


(39) 


(40) 


(41) 


and from equation (35), 


Az p = AA'e' 


leading to, 


£ 


P 

n + 1 


= £ p + AA'e' 


(42) 


(43) 


The essential step in determining the state at the next increment n + 1 has thus been reduced to the 
determination of the modified proportionality constant AA' . To do so, we write the Prandtl-Reuss 
equations, 


Ae f = Sy,n + 1 

and multiply each side by 2/3 times itself and take the square root, 


(44) 


]J 3 A8 f A8 f ~]j 3 S ij,n + 1 S >j,n + 1 ^ 


(45) 


From equation (4) and the definition of the effective stress, a = ^3/2 , equation (45) can be written as, 




(46) 


Solving equation (46) for AA and substituting into equation (38) yields. 
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E "=( l+ JSri a ^ 


(47) 


or, 


ASp 


°n + 1 

3 G 


(48) 


Then, substituting equation (40) into equation (48), we have, 


AX' = 1 - - 


3Gs 


(49) 


et 


where a is taken from the yield function, i.e., for linear hardening, 

a = Y + Hs„ 


(50) 


4. Equivalence of the Radial Return and Mendelson Methods 

It is now possible to show that the Mendelson method is exactly equivalent to the radial return 
method. First, it follows from equation (1) (and can be clearly seen by comparing equation (50) to 
equation (22)), that Mendelson’ s a is equivalent to the function k in the radial return method, which 
defines the material hardening law. That is, 


<* = k(s p ) (51) 

allowing equation (49) to be written as, 

/qe„ ) 

AX' = 1 (52) 

3 Gs et 

From the material elastoplastic constitutive equation (7), at the known increment n, we have, 

s„=2 G(e„-e*) (53) 

Substituting equation (53) into equation (9) and comparing to equation (32) yields, 

s «+i = 2(j(e„ +1 ) -2Ge' (54) 

and we see the equivalence between the role of the trial stress in the radial return method and the 
modified total strain deviator in the Mendelson method. The latter is, in fact, a strain-like trial quantity. 
This can be shown independently by considering equation (31) in the trial state, where the elastic strain 

deviator becomes the trial strain deviator and it is assumed that, in the trial condition, Az p = 0 . Equation 
(31) can then be written as, 
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(55) 


®H + 1 



+ £ 


P 

n 


T 

where the trial strain deviator is denoted as e )J+] . Rearranging equation (55) and comparing to 
equation (32), we have, 


e »+l e 


n+\ 


-rP =, 


(56) 


thus confirming that Mendelson’s modified total strain deviator is indeed a trial strain deviator in the 
sense of the radial return method. 

From equations (54), (37), and (3), 





V e ij 


\[6Gz et 


(57) 


which gives, 


V6Gl 


s »+ 1 


(58) 


and we see that Mendelson’s equivalent modified total strain takes on a role analogous to that of the norm 
of the trial stress in the radial return method. We also note that the radial return method unit normal. 


n — s «+l/ 

(57) as, 


'n+\ 


can be written in terms of Mendelson’s modified quantities using equations (54) and 


n = 



Substituting equation (59) into equation (52) yields. 



Comparing equation (18) with equation (41), we have, 



Substituting equations (58) and (60) into equation (61) gives, 


(59) 


(60) 


(61) 



( 62 ) 


NASA/TM— 2006-2 14331 


13 



Rearranging equation (62), we have, 



which is identical to the radial return method equation (17) that is employed as the consistency condition 
to determine yA t . Note that the Mendelson method equation (43) must be equivalent to the radial return 
method equation (19), which provides another route to establishing equation (63) through the Mendelson 
method equations. Thus, with the arrival at equation (63), along with the equivalence between equations 
(18) and (41) and equations (19) and (43), it has been demonstrated that the radial return and Mendelson 
methods are indeed equivalent. In fact, Mendelson’s method may be considered a total strain formulation 
of the radial return method. Conversely, of course, the radial return method may be considered a stress 
formulation of Mendelson’s total strain method. 

A final issue to consider when evaluating the equivalence between the two methods is the two 
possible alternative approaches that may be taken in the radial return method to calculate the deviatoric 
stress for increment n+ 1, step 6 in the procedure summarized in section 2. The value of s„ +1 can be 
determined from equation (14), which is derived from the yield criterion, equation (1). Alternatively, s n+1 
can be determined from the constitutive equation (53), applied at increment n+ 1, 

Sn+t = 2G(e n+1 -e£ +1 ) (64) 


as both e„ +1 and zf t+l are known at this point in the procedure. As discussed in step 5 of the radial return 
method procedure, in the case of nonlinear hardening and in the case of blended stress/strain loading, 
global iteration is needed to ensure that these two values of s n+1 are equivalent. In the context of the 
Mendelson method, it is also possible to determine the deviatoric stress components directly from the 
yield function. We first substitute equation (42) into equation (44), resulting in, 


AX' e' = s ))+ | AA, 

and from equation (39) we have, 


2G(l-AM) 

Substituting equation (66) into equation (65) yields, 

s «+i =2G(l- AA/)e' 

and substituting equation (51) into equation (67), we arrive at, 



(65) 


( 66 ) 


(67) 


( 68 ) 


which enables the calculation of the stress deviator components directly from the yield function, 
analogous to equation (14) in the radial return method. 
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To summarize the Mendelson method in a form analogous to that presented in section 2 for the radial 
return method, the process is as follows (see also table 1), 

(1) Compute elastic trial strain deviator (referred to as the modified total strain deviator by 
Mendelson), e' , from equation (32), 

p ^ — p o P 

c ^n+ 1 c n 

and compute the equivalent modified total strain, & et , from equation (37), 


8 et \l e ij e ij 


(2) Check the yield criterion, which, from equations (10) and (58), can be written as, 

/( e '> K ) = Bet = 0 

If /(e',K) <0, set s„ +1 =2Ge', e^ +1 =e^, and the increment is completed. Otherwise, 
continue. 


(3) Calculate the modified proportionality constant, AA/ , from equation (52), 

AX' = 1 - K ( £p ’” +1 ^ 


3Ge 


et 


See sections 4. 1 and 4.2 below for procedures involving linear and nonlinear hardening. 

(4) Compute the equivalent plastic strain at increment n + 1 from equation (41), 

8 P, +1 =Zp m +AX' &e t 

and the plastic strain components from equation (43) as follows, 


E«+l = £ » +AA'e' 


(5) If the hardening is nonlinear, check for global convergence by comparing the value of a„ + j 
determined from the yield function, a ;J+1 - 1 < [ z p +i j = 0 , to the value of a„ +1 determined from 
the constitutive equation (7). That is, the value of a„ +1 from the yield function is given by, 


- Y p 

a 7! + l 



(69) 


and the value of a ;j+1 from the constitutive equation (4) (utilizing equation (3)) is, 
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< -bi+l 2 [^ij,n+\ ^ij,n+\ j( e z)',«+l s ;j',H+ 1 ) 


YJ7 p o 

If cr, i+1 - g„ +1 > TOL 8 go to step 2. Note that in the case of linear hardening with loading involving 
anything other than six specified strain components, an additional convergence criterion is needed to 
ensure that the applied a„ +1 is equal to a„ +1 and a„ +i . 

(6) If necessary, calculate the deviatoric stresses from equation (68) or equation (7), 


2 /_ \ e' 

,« + i=7 k K, + i)— or s « 


= 2G \t n+x -£ p n+{ 


4.1 Linear Hardening 

In the case of linear hardening, a closed-form solution exists for determining the modified 
proportionality constant, AA' , for a given increment. This is analogous to the closed form solution for 
yA t obtained for the radial return method in section 2.1. For linear hardening within the Mendelson 
method, equations (22) and (23) are applicable. Substituting equations (22) and (41) into equation (23), 
we have. 


K ( E P,n + i) = Y + HE P,n +HA ^' E e 


Substituting equation (71) into equation (52) yields. 


AA' = 1-- 


Y + H AA' 


Solving equation (72) for AA' , we arrive at, 


1 Y + Hs 

AA' = 1 ^ 

1 + H3G 3G& ot 


4.2 Nonlinear Hardening 

In the general case of nonlinear hardening the value of AA' is obtained employing the Newton- 
Raphson algorithm. We first write equation (52) as, 

g(AA') = 1- K ( £p, ' ,+1 ) _ AA' = 0 (74) 


The procedure is then as follows (see also table 2): 
(1) Calculate for iteration k, 


o # p I p 

8 It is suggested to use a TOL value that is a small fraction of a„ + j , such as TOL = 0.0000001 ct ; j+ 
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s k+l = s k 
P> n + 1 P>n + 1 


+ (AA') 






(2) Calculate the derivative of the function g (with respect to AA' ), as follows, 
From equation (41), 


& p, n + 1 + ^ S(?? 


Expanding the function k as a lst-order Taylor series, 

K ( S „ +1 ) = K (*P» + AV £ e') = * (S,„ ) + S e , K' (e p>fl ) 

Flence, from equation (74), 


8 ( A^ 1 ' ) = G - [k ( s pn ) + AA' ' s et k' ' ( e pn ) 

■3£ e t 


- G AA ' = 0 


Therefore, 


S' = -3 K '( s /J- G 

(3) By the Newton-Raphson method. 


(AA'f +1 =(AA/)*-— 


r\k 


g 

(AA ')* 

g’ 

(AA')* 


(4) If 


g 


(AA') 


>\ k 


> TOL ” then k — ► k+ 1, go to (1) 


5. Results and Discussion 

In order to illustrate that the radial return and Mendelson Methods, as outlined above, provide 
identical results, the two methods are compared for a number of cases and the convergence behavior of 
the methods is highlighted. For a monolithic (unreinforced) elastoplastic material, we consider both pure 
strain and blended stress/strain loading conditions for both linear and exponential isotropic hardening. 
These results were generated using the FORTRAN program provided in the appendix. In addition, results 
are presented for a composite material whose elastoplastic matrix response is modeled using the radial 
return and Mendelson methods. These results were generated by implementing the two methods within 
the High-Fidelity Generalized Method of Cells (HFGMC) micromechanics model (Aboudi et al., 2003) so 
as to represent the matrix elastoplastic response. 

In all cases, the material properties (elastic modulus, Poisson’s ratio, and yield stress, respectively) 
are given by, 


E= 55.16 GPa v=0.3 7=90 MPa 


A small value for the Newton-Raphson tolerance, such as 0.0001, is suggested. 
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The material parameter, H , which corresponds to the hardening slope in the case of linear hardening (see 
eq. ( 22 )) and the initial hardening slope in the case of exponential hardening (see eq. ( 27 )) is taken as zero 
in the special case of elastic-perfectly plastic behavior. In all other cases, a value of H = 10 GPa is 
employed. The remaining material parameter required for exponential hardening, A, is varied in the 
results presented below. 

We begin by considering the linearly hardening, monolithic (unreinforced), elastoplastic material 
described above. Given the fact that classical plasticity is applicable to a volume element of material, the 
corresponding stress and strain components at this material point must be specified. Therefore, within the 
FORTRAN program, any combination of an or Sn, a 22 or 822, a33 or 833, a 2 3 or 823, an or 813, and an or 
812 may be prescribed. It should be noted that prescribing these stress and strain components within the 
computer program, often referred to as “loading” by the constitutive modeling community, is not the 
same as loading imposed on a physical boundary or surface of a structure, either in experiments or in FEA 
simulations. The first case considered involves specification of all six strain components, while all six 
stress components are unspecified and thus dictated by the elastoplastic constitutive model. Henceforth 
this loading condition will be referred to as “pure strain specification”. Note that within standard 
commercial finite element codes, this pure strain specification condition is active at the integration points 
within the elements. It is thus this type of loading condition for which the radial return method was 
originally intended. 

The normal stresses induced due to applied loading of the form 81 1 = 0 . 02 , S22 = S33 = - 0 . 01 , 823 = 813 = 
812= 0 are plotted vs. the applied Sn in figure 1 . Note that, due to the material’s isotropy, a 2 2 = a33 and a23 
= an = an = 0 . Results are shown for cases in which the loading has been applied using both 200 and 4 
increments. As Fig. 1 indicates, results from the radial return and Mendelson methods are identical for 
both the elastic-perfectly plastic (H= 0 ) and linear hardening cases irrespective of the number of 
increments employed to apply the loading. Further, the number of increments has no effect on the results 
as the solution at each of the four increments matches exactly with the corresponding 200 increment 
solution. 

Figure 2 illustrates the convergence of the 4 increment, linear hardening solution shown in figure 1 . 
Plotted is the effective stress at the first of the 4 increments (corresponding to an applied strain level of 
811 = 0 . 005 , S22 = 833 = - 0 . 0025 , 823 = S13 = S12 = 0 ) where the effective stress, a , has been calculated in 
three different ways. First, the yield function enables the calculation of a from equation ( 69 ) at each 
global iteration, m. 



( 75 ) 


Note that equation ( 75 ) can also be written from equation ( 20 ) in the radial return method where, from 
equation ( 1 ), a = y] 3/2 ||s|| . Second, the constitutive equation enables the calculation of a via 
equation ( 70 ) at each global iteration, m. 
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Figure 1. — Comparison of the Radial Return and Mendelson 
methods for a linearly hardening material (eq. (22)) under pure 
strain specification (applied strains: Sn = 0.02, 822= 833 = - 0.01, 
823= 813= 812= 0). The induced normal stresses (an and a 2 2 = 
033) are plotted vs. the applied 81 1 for different values of the 
material parameter//. Note that H= 0 corresponds to the elastic- 
perfectly plastic case. 
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Figure 2. — Comparison of the convergence behavior of the 
effective stress, a , for the Radial Return and Mendelson 
methods at an applied strain level of Sn = 0.005, s 2 2= s 33 = 

- 0.0025, s 23 = s i3 = Sn = 0 with linear hardening and H = 10 
GPa. The effective stress can be computed in three ways; from 
the material constitutive equation, from the yield function, and 
from the applied loading. All must be equal in order to ensure 
convergence. 
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Note that, again, equation (76) can be written from equation (21) in the radial return method. Finally, 

from the applied loading, the effective stress can be calculated from a = ^3/2 s“ pp s“ pp , where s app 

corresponds to the stress deviator components resulting from the applied loads. In terms of strains, this 
can be written, 




{ e ij,n+ 1) 


m I „ \m - 1 

^ij ,n + 1 


W-r-K^r 1 ' 


( 77 ) 


/ \m= 0 /_ \m /_ 

where s| „ +1 = . In the case of pure strain specification, a c ‘ pp x will correspond to a 


-c 

n + 1 


m-1 


because the strain components, and thus the strain deviator components, eu , are specified and thus do not 


change from iteration to iteration. In general, however, this is not the case. When at least one stress 
component is specified, the strain components are redistributed as a result of accumulated plastic strain 
from iteration to iteration. Equation (77) is thus independent of equation (76) in the general case. 

As shown in figure 2, for the case of pure strain specification with linear hardening, not only does 



, but also 


_r \ m l—YF \ m l—C \" !+1 

a «+i = a ;7+1 = a, ;+1 . That is, the effective stress calculated 


from the constitutive equation and the yield function are identical and do not change as a function of 
global iteration. Therefore, for this special case (pure strain specification with linear hardening), global 
iteration is not necessary and the results plotted in figure 1 can be obtained with one pass through the 
integration algorithms presented in sections 2 and 4. Note also in figure 2 that the radial return and 
Mendelson methods again yield identical results. 

In figure 3 results compare the radial return and Mendelson methods for the monolithic elastoplastic 
material with exponential hardening. As in figure 1 , the normal stresses induced due to applied loading of 
the form Sn = 0.02, 822= 833 = - 0.01, 823 = 813= 812= 0 are plotted versus the applied 81 1 for both 200 and 
4 increments. The parameter H= 10 GPa and the value of the parameter A (see eq. (27)) is varied between 
1 (which corresponds to nearly linear hardening) and 1000. Figure 3 illustrates that, once again, the radial 
return and Mendelson methods are identical and that the results are not affected by the number of 
increments employed to apply the loading. In addition, the effect of the parameter A on the exponential 
hardening results is clearly shown. Finally, note that, since all results in figure 3 employ the same value of 
the parameter H, their initial post yield slopes are identical. 

Figure 4 shows a plot of the three effective stress values vs. global iteration number for the case of 
exponential hardening with A = 250 and 4 increments at the first increment. As was the case in figure 2, 


figure 4 shows that 



, which will always be the case in pure strain specification. 


Now, however, the effective stress values calculated from the yield function and from the constitutive 
equation do not correspond until convergence has been achieved. This is because of the nonlinear 

function employed to describe the material’s hardening behavior. Note that convergence is rapid; a„ +1 

and a n+1 are within 1 percent of each other after two iterations. 

Next, consider a loading condition involving specification of both stress and strain components on the 
elastoplastic material. This loading condition will be referred to as “blended specification”. The loading 
will take the form of a specified normal strain component in one direction with a stress-free condition 
specified for all other components. That is, we apply, sn = 0.02, G 22 = G 33 = c ?23 = G 13 = G 12 = 0, which 
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Figure 3. — Comparison of the Radial Return and Mendelson 
methods for a material represented by an exponentially 
hardening material (eq. (27)) under pure strain specification 
(applied strains: Sn = 0.02, £22= 833 = - 0.01, 823= £13 = 812= 0). 
The induced normal stresses (an and O 22 = a33) are plotted vs. 
the applied 81 1 for different values of the material parameter A, 
with H= 10 GPa. 
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Figure 4. — Comparison of the convergence behavior of the 
effective stress, a , for the Radial Return and Mendelson 
methods at an applied strain level of Sn = 0.005, 822= s 33 = 

- 0.0025, £23 = 813 = 812 = 0 with exponential hardening and A = 
250. The effective stress can be computed in three ways; from 
the material constitutive equation, from the yield function, and 
from the applied loading. All must be equal in order to ensure 
convergence. 
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simulates the state experienced by the material in a standard strain controlled tensile test. Figure 5 plots 
the stress-strain curves in the applied strain direction (i.e., an vs. Sn), along with the induced normal 
strain in the unapplied direction (i.e., an vs. 822 = 833) for a linear hardening elastoplastic material where 
the loading has been applied using 200 and 4 increments. The radial return and Mendelson methods are 
compared. As in figures 1 and 3, figure 5 shows that the radial return and Mendelson methods provide 
identical results in both the elastic-perfectly plastic (H= 0) and linear hardening cases. The number of 
increments employed also has no effect on the results. 

Figure 6 provides a plot of the global convergence behavior of the effective stress (again calculated in 
three different ways) for the first increment of the 4 increment case plotted in figure 5. As was shown in 
figures 2 and 4, the convergence of the radial return and Mendelson methods are again identical. Further, 
as was the case in figure 2 , since the material hardening is linear, the effective stresses calculated from the 
yield function and from the constitutive equation are identical. Flowever, due to the blended loading 
specification, the effective stress calculated from the applied loading is independent. This is because, in 
the case of blended specification, the strain deviator components change from iteration to iteration and 
thus equations (77) and (76) are independent. Recall that under pure strain specification, the strain 
deviator components are known exactly for a given load increment and do not change as a function of 

\m / si \m— 1 

®n P +\ = ®n + 1 • The key point illustrated in figure 6 is that, 

because the yield function and constitutive equation effective stress values are identical, but not 
converged, these two effective stress values cannot be used as a measure of convergence in the case of 
linear hardening with loading other than pure strain control. Were such a measure used, both the radial 
return and Mendelson methods would converge immediately to an incorrect solution that violates the 
specified loading condition (e.g., the stress components prescribed as zero would be non-zero). 



811 , 822 

Figure 5. — Comparison of the Radial Return and Mendelson 
methods for a material represented by a linearly hardening 
material (eq. ( 22 )) under blended specification (applied load: 
8 n = 0.02, a 22 = a 33 = 023 = an = an = 0). The normal strains 
(sn and S 22 = S33) are plotted vs. an for different values of the 
material parameter//. Note that H= 0 corresponds to the elastic- 
perfectly plastic case. 
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Figure 6 . — Comparison of the convergence behavior of the 
effective stress, a , for the Radial Return and Mendelson 
methods at an applied strain level of Sn = 0.005, a 22 = <733 = 033 = 

( 7 i 3 = a 12 = 0 with linear hardening and H= 10 GPa. The 
effective stress can be computed in three ways; from the material 
constitutive equation, from the yield function, and from the 
applied loading. All must be equal in order to ensure 
convergence. 

The radial return and Mendelson methods are compared in figure 7 in the case of blended 
specification for an exponential hardening material with a varying value of the parameter A. The specified 
loading is identical to that employed in figure 5, and 200 and 4 increments have again been used. As 
observed previously, the radial return and Mendelson methods provide identical results and the number of 
increments employed is of no consequence. In figure 8 , the convergence behavior of the effective stress is 
shown for the first increment of the 4 increment case from figure 7 with A = 250. Once again, the radial 
return and Mendelson methods are identical. As was the case in figure 6 , the effective stress calculated 
from the applied loading is independent of that calculated from the constitutive equation due to the 
blended specification (see eqs. (76) and (77)). Unlike figure 6 , however, figure 8 indicates that, due to the 
nonlinear hardening behavior of the elastoplastic material, the effective stress calculated from the 
constitutive equation and the yield function are no longer coincident. Thus, in this case, these two 
effective stress values could be used as a measure of global convergence. 

As a final illustration of the radial return and Mendelson methods, both implementations have been 
incorporated into a micromechanics model to enable the prediction of the response of a continuous 
composite material with elastoplastic phases. The High-Fidelity Generalized Method of Cells (HFGMC) 
micromechanics model (Aboudi et. al, 2003) has been employed for this purpose. This semi-analytical 
(non-FEA) model discretizes the composite repeating unit cell geometry into an arbitrary number of 
subcells, each of which may contain a distinct, arbitrary material. HFGMC localizes an arbitrary, 
admissible combination of stress and strain components on a composite material to determine the point- 
wise stresses and strains throughout the fiber and matrix constituents. These point-wise stress and strain 
components can thus be passed to a local constitutive model (such as classical plasticity) in order to 
simulate the local plastic response at each point. The local plastic strains are then homogenized by 
HFGMC to determine composite-level inelastic strains. Note that, in the linear elastic case, the HFGMC 
localization/homogenization method provides the effective thermo-elastic properties of the composite 
materials from those of the constituents. 
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811 , 822 

Figure 7 . — Comparison of the Radial Return and Mendelson 
methods for a material represented by a exponentially hardening 
material (eq. ( 27 )) under blended specification (applied load: Sn = 
0 . 02 , 022= O33 = 023= 013= 012= 0 ). The normal strains (sn and 822 
= 833) are plotted vs. On for different values of the material 
parameter A, with 77 = 10 GPa. 



Figure 8. — Comparison of the convergence behavior of the effective 
stress, o , for the Radial Return and Mendelson methods at an 
applied strain level of Sn = 0 . 005 , 022= c?33 = 023= 013= 012= 0 
exponential hardening and A = 250 . The effective stress can be 
computed in three ways; from the material constitutive equation, 
from the yield function, and from the applied loading. All must be 
equal in order to ensure convergence. 
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The HFGMC repeating unit cell employed herein is shown in figure 9, where the fiber and matrix 
constituents are depicted in different colors. As shown, in the plane of the fiber’s cross-section, the 
repeating unit cell is discretized into 8 by 8 subcells, whereas, in the out of plane direction, the material is 
considered to be infinitely long. The matrix material’s properties E, v, and Y are taken as those of the 
elastoplastic material considered above. The fiber material is taken to be elastic with properties given by, 

E = 413.7 GPa v=0.2 

We consider blended stress/strain component specification on a 35 percent fiber volume fraction 
composite, simulating a transverse uniaxial strain controlled tensile test such that s 2 2 = 0.01, an = a 33 = 
a 23 = c?i 3 = an = 0. Results are shown in figure 10, where the composite stress-strain curves in the applied 
direction (i.e., a 22 vs. s 22 ) are plotted, along with the normal strains induced in the composite (i.e., a 22 vs. 
Si i and a 22 vs. s 33 ). Three types of elastoplastic matrix materials are considered: elastic-perfectly plastic, 
linear hardening with H = 10 GPa, and exponential hardening with H = 10 GPa and A = 100. Cases 
involving 1 00 increments and 4 increments are considered, and the radial return and Mendelson methods 
are once again compared. First, it is clear from figure 10 that the stiff continuous fiber in thexi-direction 
restrains the composite in this direction, resulting in an induced Sn that is significantly smaller than the 
induced s 33 . Second, we see that, as in all other cases, the radial return and Mendelson methods produce 
identical results. Unlike previous results, however, the number of increments employed does have an 
effect on the predicted composite response, with the results for 4 increments slightly underpredicting the 
composite stress-strain curves. It was determined (but not shown) that for the present case, approximately 
10 increments were required to achieve good convergence of the composite level stresses and strains. 



Figure 9. — HFGMC micromechanics 
theory repeating unit cell employed to 
simulate the transverse behavior of a 
composite material. 
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811 , 822 , 833 

Figure 10. — Comparison of the Radial Return and Mendelson 
methods utilized to simulate the matrix behavior within a 
35 percent volume fraction composite material under blended 
specification (applied load: s 22 = 0.02, an = a 33 = cr 23 = a i3 = a 12 = 
0; transverse loading). The normal strains (sn, s 22 , and s 33 ) are 
plotted vs. a 22 for linear hardening, exponential hardening, and 
elastic-perfectly plastic matrix behavior. 


Figure 1 1 illustrates sample convergence behavior for the composite simulations. As in figures 2, 4, 6, 
and 8, figure 1 1 illustrates the effective stress convergence for the first increment in the case that 
employed 4 increments. Figure 1 1 is further specific to the exponential hardening case from figure 10. In 
addition, because the local behavior is different at each location within the composite matrix, figure 1 1 is 
applicable to one particular matrix location; that specified by the red “x” in figure 9. Figure 1 1 shows that 
the radial return and Mendelson methods again provide identical results. Further, as was the case for the 
monolithic (unreinforced) elastoplastic material with blended specification, the effective stress calculated 
from the constitutive equation and the load/localization procedure are unique. In fact, in the case of a 
composite, because the stresses and strains are localized from the applied composite-level stresses and 
strains, the local convergence will behave in this way regardless of the form of the specified composite 
level loading. That is, the distribution of the strain components (and thus the strain deviator components, 
e,j) at a given point are always free to change from iteration to iteration regardless of the form of the 
specified loading on the composite. This is because this local strain component distribution is dictated by 
the micromechanics localization, which is influenced by the state at every other point in the composite, 
each of which is at a different point along its yield function. It is also noteworthy in figure 1 1 that the 
convergence is considerably more gradual compared to the monolithic elastoplastic material convergence 
illustrated in figures 2, 4, 6, and 8. 
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Figure 1 1 . — Comparison of the global convergence behavior of the 
effective stress, a , for the Radial Return and Mendelson methods 
within a 35 percent volume fraction composite material (at the 
location indicated in figure 9) at an applied global (composite) 
strain level of 822 = 0.0025, (Jn = a 33 = a 23 = CTi 3 = CTi 2 = 0 with 
exponential hardening and A = 100. The effective stress can be 
computed in three ways; from the material constitutive equation, 
from the yield function, and from the applied 
loading/micromechanics localization. All must be equal in order to 
ensure convergence. 


6. Conclusions 

The Mendelson (1968, 1983) and radial return (Simo and Taylor, 1985; Simo and Flughes, 1998) 
methods for the integration of the classical plasticity are equivalent. Key to recognizing this fact is the 
identification of the equivalence between the role of the trial stress in the radial return method and that of 
Mendelson’ s modified total strain deviator, with the latter being a strain-like trial quantity. The 
Mendelson method may thus be thought of as a total strain formulation of the radial return method. The 
algorithmic implementation of the methods is very similar, and the results of both methods are identical. 
As such, there appears to be no advantage computationally or theoretically to one method or the other. 
The radial return method may have a conceptual advantage in that returning to the yield surface from the 
trial stress state can be visualized. The steps involved in the two methods are summarized in tables 1 and 2. 

Results for the radial return and Mendelson methods have been compared for the cases of linear and 
nonlinear isotropic hardening of a monolithic (unreinforced) elastoplastic material. The 
articulation/implementation of explicit nonlinear hardening within the Mendelson method presented 
herein is similar to that provided by Simo and Taylor (1985) for the radial return method. In addition to 
illustrating the numerical equivalence of the two integration methods (in all cases), the presented results 
focus on the effects of the type of loading (pure strain specification or blended stress/strain specification), 
number of increments employed to apply the loading, and the convergence behavior of the methods as 
local iteration is performed to ensure that the plastic state remains on the yield surface. The following 
conclusions, applicable to both the radial return and Mendelson methods, can be drawn from the 
monolithic material results presented: 
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(1) The results are independent of the number of increments used to apply the load. 

(2) Local iteration is not necessary in the case of pure strain specification with linear isotropic 
hardening. 

(3) Local iteration is necessary in the cases of blended specification or nonlinear isotropic hardening. 

(4) Blended specification with linear isotropic hardening represents a special case in that the effective 
stress value calculated from the constitutive equation and the yield function are identical for a 
given iteration, though not necessarily converged. These two values thus cannot be used as a 
convergence criterion in this case. In this case, the effective stress value calculated from the 
applied loading may be compared with one of the other two effective stress values to gauge 
convergence. However, employing this applied loading effective stress value in all cases to 
determine convergence may result in the performance of an extra, unnecessary iteration. 

(5) Pure strain specification behaves differently than blended specification because the distribution 
(i.e., proportion) of the strain components is prescribed and thus known a priori. The plasticity 
equations then need only separate each strain component into elastic and inelastic parts. In 
contrast, under blended specification, the plasticity equations can find a strain component 
distribution corresponding to a plastic state that is on the yield surface, but does not satisfy the 
specified blended loading. 

Results for a composite material whose elastoplastic matrix material was modeled using the radial 
return and Mendelson classical plasticity methods were also presented. The composite analysis employed 
the non-FEA High-Fidelity Generalized Method of Cells (HFGMC) micromechanics model (Aboudi et 
al., 2003) developed by the authors. Results using both integration methods were once again identical, 
however, in contrast to the monolithic material results, the number of increments employed to apply the 
loading now had an effect. Finally, because the local stresses and strains have been localized from the 
applied composite-level stresses and strains, the convergence behavior is similar to that of the monolithic 
case with blended specification and nonlinear hardening. Iteration is required regardless of the nature of 
the loading applied on the composite and the type of hardening. This fact is somewhat intuitive because 
the in situ stresses within the micromechanics representation of the composite are multi-axial and 
progress in a non-proportional way irrespective of the nature of the global composite-level loading. Thus, 
on the local level, the micromechanics problem is similar to the blended specification monolithic material 
case where the correct proportional distribution of the strain components is not known a priori. 
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Appendix 


Radial Return and Mendelson Methods FORTRAN Program 


PROGRAM MAIN 


-- Integrates the classical plasticity equations for arbitrary monotonic loading using the Radial Return and 
Mendelson methods for linear and exponential isotropic hardening materials 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

CHARACTER* 3 : : HARDENING 
INTEGER :: LOP (6) 

DOUBLE PRECISION :: APPLYF(6) 


! — Properties 
E = 55160. 

XNU =0.3 
YIELD = 90. 
HARD = 10000. 

! HARD = 0. 

A = 1000. 

! A = 0.0000001 


! — 6 Components of applied loading 
! — LOP = 1 > strain; LOP = 2 > stress 


LOP ( 1 ) = 1 


LOP (2) = 2 


LOP (3) = 2 


LOP (4) = 2 


LOP (5) = 2 


LOP (6) = 2 


APPLYF (1) = 

0.02 

APPLYF (2) = 

0 . 

APPLYF (3) = 

0 . 

APPLYF (4) = 

0 . 

APPLYF (5) = 

0 . 

APPLYF (6) = 

0 . 


! -- Number of increments 
NINC = 4 

! HARDENING = 'LIN' 

HARDENING = 'NON' 

TOL_GLOBAL = 0.0000001 
TOL_NEWTON = 0.0001 

CALL RADIAL_RETURN (E, XNU, YIELD, HARD, A, LOP, APPLYF, TOL_GLOBAL, TOL_NEWTON, NINC, HARDENING) 
CALL MENDELSON (E, XNU, YIELD, HARD, A, LOP, APPLYF, TOL_GLOBAL, TOL_NEWTON, NINC, HARDENING) 

END 


SUBROUTINE RADIAL_RETURN (E, XNU, YIELD, HARD, A, LOP, APPLYF, TOL_GLOBAL, TOL_NEWTON, NINC, HARDENING) 


— Integrates the classical plasticity equations for arbitrary monotonic loading 
using the Radial Return method 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

CHARACTER* 3 : : HARDENING 
LOGICAL : : CONVERGED 
INTEGER :: LOP (6) 

DOUBLE PRECISION :: APPLYF (6), APPLY (6) 


OPEN (22, FILE= ' radial_return . plot ' , ST ATUS= ' UNKNOWN ' ) 
OPEN (23, FILE= ' radial_return . conv ' , ST ATUS= ' UNKNOWN ' ) 
OPEN (24, FILE= ' radial_return . out ' , ST ATUS= ' UNKNOWN ' ) 

WRITE (24, *) 'RADIAL RETURN METHOD OUTPUT' 

WRITE (24, *) ' ' 
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WRITE (24, 9700) E, XNU 
IF (HARDENING .EQ. 'LIN') THEN 

WRITE (24, *) 'LINEAR HARDENING SELECTED' 

WRITE (24, 9710) YIELD, HARD 
ELSE 

WRITE (24, *) 'EXPONENTIAL HARDENING SELECTED' 
WRITE (24, 9720) YIELD, HARD, A 
ENDIF 

WRITE (24, 9725) LOP 
WRITE (24, 9730) APPLYF 
WRITE (24, 9735) NINC 

WRITE (24, 9740) TOL_GLOBAL, TOL_NEWTON 

C44 = E/ (2.* (1. + XNU) ) 

BULK = E/(3.*(l. - 2*XNU) ) 

SQ2 = DSQRT (2 . DO) 

SQ3 = DSQRT (3. DO) 


EPS11 = 0 
EPS22 = 0 
EPS33 = 0 
EPS23 = 0 
EPS13 = 0 
EPS12 = 0 

E11P = 0 
E22P = 0 
E33P = 0 
E23P = 0 
E13P = 0 
E12P = 0 

! -- Incremental loading loop 

DO 1000 WHILE (INC . LT . NINC) 

INC = INC + 1 

WRITE (23, *) 'INC = ' , INC 

! -- Increment applied loading 
DO I = 1, 6 

APPLY (I) = APPLY (I) + APPLYF (I) /NINC 
END DO 


I GLOB = 0 

CONVERGED = .FALSE. 

! — Global iteration loop 

DO 2000 WHILE (.NOT. CONVERGED) 

I GLOB = I GLOB + 1 

! -- Determine stresses and strains from applied loading 

CALL LOAD ( E , XNU, C44, LOP, APPLY, EPS11, EPS22 , EPS33, EPS23, EPS13, EPS12 , & 
EllP, E22P, E33P, E23P, E13P, E12P, SEFF_LOAD) 


! — Computation of the strain deviators [ EDij ] (This is at n+1) 

EPSMEAN= (EPS11 + EPS22 + EPS33) / 3 

ED11 = EPS11 - EPSMEAN 
ED22 = EPS22 - EPSMEAN 
ED33 = EPS33 - EPSMEAN 
ED23 = EPS23 
ED13 = EPS13 
ED12 = EPS12 

PRESS = BULK* EPSMEAN* 3. 

! — Computation of the trial stress deviators (STij) 

ST11 = 2*C44* (ED11 - EllP) 

ST22 = 2*C44* (ED22 - E22P) 

ST33 = 2*C44* (ED33 - E33P) 

ST23 = 2*C44* (ED23 - E23P) 

ST13 = 2*C44* (ED13 - E13P) 

ST12 = 2*C44* (ED12 - E12P) 

! -- Computation of the unit normal vector (XNij) 

XNORM = DSQRT (ST11**2 + ST22**2 + ST33**2 + 2*(ST23)**2 + 2*(ST13)**2 + 2*(ST12)**2) 

! — Elastic case (inside initial yield surface) 

IF ( XNORM* SQ3/SQ2 . LT . YIELD) THEN 
SIG11C = ST11 + PRESS 

SIG22C = ST22 + PRESS 

SIG33C = ST33 + PRESS 

SIG23C = ST23 

SIG13C = ST13 

SIG12C = ST12 

GOTO 999 
ENDIF 
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XN11 = ST11 /XNORM 
XN22 = ST22/XNORM 
XN33 = ST33/XNORM 
XN23 = ST23/XNORM 
XN13 = ST13/XNORM 
XN12 = ST12/XNORM 


! — Calculate GAMMADT 


! -- Linear hardening 

IF (HARDENING .EQ. 'LIN') THEN 

GAMMADT = (XNORM - (SQ2/SQ3) * (YIELD + HARD*EPQ) ) / (2*C44 + 2.*HARD/3.) 


! — Exponential hardening 
ELSE 

INEWT = 0 
GAMMADT = 0 . DO 

10 CONTINUE 

EPQ1 = EPQ + (SQ2/SQ3) * GAMMADT 

XKAPPA = YIELD - (HARD/A) *( DEXP ( -A*EPQ1 ) - 1) 

DXKAPPA = HARD*DEXP (-A*EPQ1) 

G = - (SQ2/SQ3) *XKAPPA - 2 . D0*GAMMADT*DXKAPPA/3 . 0 + XNORM - 2 . D0*C44*GAMMADT 
DG = -2. DO* DXKAPPA/ 3. DO - 2.D0*C44 

GAMMADT = GAMMADT - G/DG 

INEWT = INEWT + 1 

IF (DABS (G) .GT. TOL_NEWTON) GOTO 10 
ENDIF 


EPQ = EPQ + (SQ2/SQ3) * GAMMADT 


E11P = E11P + 
E22P = E22P + 
E33P = E33P + 
E23P = E23P + 
E13P = E13P + 
E12P = E12P + 


GAMMADT *XN 11 
GAMMADT *XN2 2 
GAMMADT *XN3 3 
GAMMADT *XN2 3 
GAMMADT *XN 13 
GAMMADT *XN 12 


IF (HARDENING .EQ. 'LIN') THEN 
XKAPPA = YIELD + HARD* EPQ 
ELSE 

XKAPPA = YIELD - (HARD/A) * (DEXP (-A*EPQ) - 1) 

ENDIF 

! -- Calculate effective stress from constitutive eqn (= XKAPPA for convergence) 

SEFF = 2*C44* (SQ3/SQ2) *DSQRT ( (ED11 - E11P) **2 + (ED22 - E22P) **2 + (ED33 - E33P) **2 + 
2* (ED23 - E23P) **2 + 2* (ED13 - E13P) **2 + 2* (ED12 - E12P) 

WRITE (23, 9100) IGLOB, SEFF, XKAPPA, SEFF_LOAD, EPQ 


! -- Global convergence check 
CONVERGED = .TRUE. 

IF (DABS ((SEFF - XKAPPA) /SEFF) . GT . TOL_GLOBAL) CONVERGED = .FALSE. 


! — Additional check based on applied loading (needed for blended loading with linear hardening) 
IF (DABS ( (SEFF_LOAD - XKAPPA) /SEFF_LOAD) . GT . TOL_GLOBAL) CONVERGED = .FALSE. 


2000 CONTINUE 


Sll = (SQ2/SQ3) *XKAPPA*XN11 
S22 = (SQ2/SQ3) *XKAPPA*XN22 
S33 = (SQ2/SQ3) *XKAPPA*XN33 
SIG23 = (SQ2/SQ3) *XKAPPA*XN23 
SIG13 = (SQ2/SQ3) *XKAPPA*XN13 
SIG12 = (SQ2/SQ3) *XKAPPA*XN12 


SIGll = 

Sll + PRESS 



SIG22 = 

S22 + PRESS 



SIG33 = 

S33 + PRESS 



SIG11C = 

2*C44* (ED11 - 

E11P) + 

PRESS 

SIG22C = 

2*C44* (ED22 - 

E22P) + 

PRESS 

SIG33C = 

2*C44* (ED33 - 

E33P) + 

PRESS 

SIG23C = 

2*C44* (ED23 - 

E23P) 


SIG13C = 

2*C44* (ED13 - 

E13P) 


SIG12C = 

2*C44* (ED12 - 

E12P) 


WRITE (22 

, 9000) EPS11 , 

SIG11C, 

EPS22. 


WRITE (24, 
WRITE (24, 
WRITE (24, 
WRITE (24, 


9200) INC, IGLOB 

9300) EPS11 , EPS22 , EPS33, EPS23, EPS13, EPS12 
9400) E11P, E22P, E33P, E23P, E13P, E12P 
9500) SIGH, SIG22, SIG33, SIG23, SIG13, SIG12 
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WRITE (24, 9600) SIG11C, SIG22C, SIG33C, SIG23C, SIG13C, SIG12C 


1000 CONTINUE 

CLOSE (22) 
CLOSE (23) 
CLOSE (24) 

RETURN 


9000 FORMAT (4 (E12 . 5, 2X) ) 

9100 FORMAT (13, 2X, 4 (E12.5,2X) ) 

9200 FORMAT (//,' INCREMENT # ', 14, 6X, ' # GLOBAL ITERATIONS ’,14) 

9300 FORMAT (' TOTAL STRAIN COMPONENTS 3X, 6 (E12 . 5, IX) ) 

9400 FORMAT (' PLASTIC STRAIN COMPONENTS 3X, 6 (E12 . 5, IX) ) 

9500 FORMAT (' STRESS COMPONENTS FROM YIELD FUNCTION: ' , /, 3X, 6 (E12 . 5, IX) ) 

9600 FORMAT (' STRESS COMPONENTS FROM CONSTITUTIVE EQN: ' , /, 3X, 6 (E12 . 5, IX) ) 

9700 FORMAT (//,' ELASTIC MATERIAL PARAMETERS 3X, ' ELASTIC MODULUS ', E12 . 5, /, 3X, ' POISSON RATIO ’,F7.5,/) 

9710 FORMAT ( / , ' PLASTIC MATERIAL PARAMETERS 3X, ' YIELD STRESS ', E12 . 5 ,/, 3X, ' HARDENING SLOPE ’,E12.5) 

9720 FORMAT (/,' PLASTIC MATERIAL PARAMETERS 3X, ' YIELD STRESS ' , E12 . 5, /, 3X, ' HARDENING SLOPE ' , E12 . 5 , / , 3X, & 

'EXP. PARAMETER A ’,E12.5) 

9725 FORMAT (/,' APPLIED LOAD (1 = STRAIN; 2 = STRESS) : ' , 6 (112, IX) ) 

9730 FORMAT (5X, 6 (E12.5, IX) ) 

9735 FORMAT (/,' NUMBER OF INCREMENTS ',14) 

9740 FORMAT ( / , ' TOLERANCE FOR GLOBAL ITERATIONS: ',E12.5,/,' TOLERANCE FOR NEWTON-RAPHSON ITERATIONS: ’,E12.5) 

END 


SUBROUTINE MENDELSON (E, XNU, YIELD, HARD, A, LOP, APPLYF, TOL_GLOBAL, TOL_NEWTON, NINC, HARDENING) 


— Integrates the classical plasticity equations for arbitrary monotonic loading 
using the Mendelson method 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

CHARACTER* 3 : : HARDENING 
LOGICAL : : CONVERGED 

INTEGER :: LOP (6) 

DOUBLE PRECISION :: APPLYF (6), APPLY (6) 


OPEN (19, FILE= 'mendelson. plot ' , ST ATUS= ' UNKNOWN ' ) 
OPEN (20, FILE= ' mendelson . conv ' , STATUS= ' UNKNOWN ' ) 
OPEN (21, FILE= 'mendelson. out ' , ST ATUS= ' UNKNOWN ' ) 

WRITE (21, *) 'MENDELSON METHOD OUTPUT' 

WRITE (21, *) ' ' 


WRITE (21, 9700) E, XNU 
IF (HARDENING .EQ. 'LIN') THEN 

WRITE (21, *) 'LINEAR HARDENING SELECTED' 

WRITE (21, 9710) YIELD, HARD 
ELSE 

WRITE (21, *) 'EXPONENTIAL HARDENING SELECTED' 
WRITE (21, 9720) YIELD, HARD, A 
ENDIF 


WRITE (21, 

9725) 

LOP 

WRITE (21, 

9730) 

APPLYF 

WRITE (21, 

9735) 

NINC 

WRITE (21, 

9740) 

TOL_GLOBAL, 

C44 = E/ ( 

2 . * ( 1 . 

+ XNU) ) 

BULK = E/ (3. * (1 

. - 2*XNU) ) 


TOL NEWTON 


SQ2 = DSQRT (2 . DO) 
SQ3 = DSQRT (3. DO) 


EPS11 = 0 
EPS22 = 0 
EPS33 = 0 
EPS23 = 0 
EPS13 = 0 
EPS12 = 0 

E11P = 0 
E22P = 0 
E33P = 0 
E23P = 0 
E13P = 0 
E12P = 0 
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! -- Incremental loading loop 

DO 1000 WHILE (INC . LT . NINC) 

INC = INC + 1 

WRITE (20, *) 'INC = ' , INC 

! — Increment applied loading 
DO I = 1, 6 

APPLY (I) = APPLY (I) + APPLYF(I) /NINC 
END DO 


I GLOB = 0 

CONVERGED = .FALSE. 

! — Global iteration loop 

DO 2000 WHILE (.NOT. CONVERGED) 

I GLOB = I GLOB + 1 


! -- Determine stresses and strains from applied loading 

CALL LOAD ( E , XNU, C44, LOP, APPLY, EPS11, EPS22 , EPS33, EPS23, EPS13, EPS12 , & 
EllP, E22P, E33P, E23P, E13P, E12P, SEFF_LOAD) 

! -- Computation of the strain deviators [ EDij ] (This is at n+1) 

EPSMEAN= (EPS11 + EPS22 + EPS33) / 3 

ED11 = EPS11 - EPSMEAN 
ED22 = EPS22 - EPSMEAN 
ED33 = EPS33 - EPSMEAN 
ED23 = EPS23 
ED13 = EPS13 
ED12 = EPS12 

PRESS = BULK*EPSMEAN*3 . 


! -- Computation of the modified total strain deviators (BEPTij) (= Trial strain deviators) 


BEPT11 

= ED11 

- EllP 

BEPT22 

= ED22 

- E22P 

BEPT33 

= ED33 

- E33P 

BEPT23 

= ED23 

- E23P 

BEPT13 

= ED13 

- E13P 

BEPT12 

= ED12 

- E12P 


! — Computation of the equivalent modified total strain (bepet) 

BEPET = (SQ2/SQ3) *DSQRT (BEPT11**2 + BEPT22**2 + BEPT33**2 + 2* (BEPT23) **2 + & 
2* (BEPT13) **2 + 2* (BEPT12) **2) 

! — Elastic case (inside initial yield surface) 

IF (BEPET .LE. YIELD/ ( 3 . *C44 ) ) THEN 
SIG11C = 2*C44*BEPT11 + PRESS 

SIG22C = 2*C44*BEPT22 + PRESS 

SIG33C = 2*C44*BEPT33 + PRESS 

SIG23C = 2*C44*BEPT23 

SIG13C = 2*C44*BEPT13 

SIG12C = 2*C44*BEPT12 

GOTO 999 
ENDIF 


! — Calculate DLAM 


! — Linear hardening 

IF (HARDENING .EQ. 'LIN') THEN 

DLAM = (3*C44 - (YIELD + HARD*EPQ) /BEPET) / (3*C44 + HARD) 


! -- Exponential hardening 
ELSE 

INEWT = 0 
DLAM = 0.D0 

10 CONTINUE 

EPQ1 = EPQ + DLAM*BEPET 

XKAPPA = YIELD - (HARD/A) * (DEXP (-A*EPQ1) - 1) 

DXKAPPA = HARD*DEXP (-A*EPQ1) 

G = C44 - (XKAPPA + DLAM*BEPET*DXKAPPA) / (3 . *BEPET) - DLAM*C44 
DG = -DXKAPPA/ 3 . - C44 

DLAM = DLAM - G/DG 

INEWT = INEWT + 1 

IF (DABS (G) .GT. TOL_NEWTON) GOTO 10 
ENDIF 


EPQ = EPQ + DLAM*BEPET 
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E11P = E11P + DLAM*BEPT11 

E22P = E22P + DLAM*BEPT22 

E33P = E33P + DLAM*BEPT33 

E23P = E23P + DLAM*BEPT23 

E13P = E13P + DLAM*BEPT13 

E12P = E12P + DLAM*BEPT12 

IF (HARDENING .EQ. 'LIN') THEN 
XKAPPA = YIELD + HARD*EPQ 
ELSE 

XKAPPA = YIELD - (HARD/A) * (DEXP (-A*EPQ) - 1) 

ENDIF 

! -- Calculate effective stress from constitutive eqn (= XKAPPA for convergence) 

SEFF = 2*C44* (SQ3/SQ2) *DSQRT ( (ED11 - E11P) **2 + (ED22 - E22P) **2 + (ED33 - E33P) **2 + & 

2* (ED23 - E23P) **2 + 2* (ED13 - E13P) **2 + 2* (ED12 - E12P)**2) 

WRITE (20, 9100) IGLOB, SEFF, XKAPPA, SEFF_LOAD, EPQ 

! -- Global convergence check 
CONVERGED = .TRUE. 

IF (DABS ((SEFF - XKAPPA) /SEFF) . GT . TOL_GLOBAL) CONVERGED = .FALSE. 

! — Additional check based on applied loading (needed for blended loading with linear hardening) 

IF (DABS ( (SEFF_LOAD - XKAPPA) /SEFF_LOAD) . GT . TOL_GLOBAL) CONVERGED = .FALSE. 

2000 CONTINUE 


SIGH = 2*XKAPPA*BEPT11/ (3. *BEPET) + PRESS 
SIG22 = 2*XKAPPA*BEPT22/ (3. *BEPET) + PRESS 
SIG33 = 2*XKAPPA*BEPT33/ (3. *BEPET) + PRESS 
SIG23 = 2*XKAPPA*BEPT23/ (3. *BEPET) 

SIG13 = 2*XKAPPA*BEPT13/ (3. *BEPET) 

SIG12 = 2*XKAPPA*BEPT12/ (3. *BEPET) 

SIG11C = 2*C44* (ED11 - E11P) + PRESS 

SIG22C = 2*C44* (ED22 - E22P) + PRESS 

SIG33C = 2*C44* (ED33 - E33P) + PRESS 

SIG23C = 2*C44* (ED23 - E23P) 

SIG13C = 2*C44* (ED13 - E13P) 

SIG12C = 2*C44* (ED12 - E12P) 


999 WRITE (19, 9000) EPS11, SIG11C, EPS22, EPS33 
WRITE (21, 9200) INC, IGLOB 

WRITE (21, 9300) EPS11, EPS22 , EPS33, EPS23, EPS13, EPS12 
WRITE (21, 9400) E11P, E22P, E33P, E23P, E13P, E12P 
WRITE (21, 9500) SIGll, SIG22, SIG33, SIG23, SIG13, SIG12 
WRITE (21, 9600) SIG11C, SIG22C, SIG33C, SIG23C, SIG13C, SIG12C 

1000 CONTINUE 

CLOSE (19) 

CLOSE (20) 

CLOSE (21) 

RETURN 


9000 FORMAT (4 (E12 . 5, 2X) ) 

9100 FORMAT (13, 2X, 4 (E12.5,2X) ) 

9200 FORMAT (//,' INCREMENT # ', 14, 6X, ' # GLOBAL ITERATIONS ',14) 

9300 FORMAT (' TOTAL STRAIN COMPONENTS :',/, 3X, 6 (E12 . 5, IX) ) 

9400 FORMAT (' PLASTIC STRAIN COMPONENTS :',/, 3X, 6 (E12 . 5 , IX) ) 

9500 FORMAT (' STRESS COMPONENTS FROM YIELD FUNCTION: ' , /, 3X, 6 (E12 . 5, IX) ) 

9600 FORMAT (' STRESS COMPONENTS FROM CONSTITUTIVE EQN: ' , /, 3X, 6 (E12 . 5, IX) ) 

9700 FORMAT ( / / , ' ELASTIC MATERIAL PARAMETERS :',/, 3X, ' ELASTIC MODULUS ', E12 . 5, /, 3X, ' POISSON RATIO ',F7.5,/) 

9710 FORMAT (/,' PLASTIC MATERIAL PARAMETERS :',/, 3X, ' YIELD STRESS ', E12 . 5, /, 3X, ' HARDENING SLOPE ',E12.5) 

9720 FORMAT ( / , ' PLASTIC MATERIAL PARAMETERS :',/, 3X, ' YIELD STRESS ', E12 . 5, /, 3X, ' HARDENING SLOPE ' , E12 . 5, /, 3X, & 

'EXP. PARAMETER A ',E12.5) 

9725 FORMAT (/,' APPLIED LOAD (1 = STRAIN; 2 = STRESS) : ' , /, 6 (112, IX) ) 

9730 FORMAT (5X, 6 (E12.5, IX) ) 

9735 FORMAT (/,' NUMBER OF INCREMENTS ',14) 

9740 FORMAT ( / , ' TOLERANCE FOR GLOBAL ITERATIONS: ',E12.5,/,' TOLERANCE FOR NEWTON-RAPHSON ITERATIONS: ',E12.5) 

END 



SUBROUTINE LOAD ( E , XNU, C44, LOP, APPLY, EPS11, EPS22 , EPS33, EPS23, EPS13, EPS12 , & 
EllP, E22P, E33P, E23P, E13P, E12P, SEFF_LOAD) 


-- Given some combination of known S and E, solves for unknown components from S = C* (E - EP) 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

INTEGER :: LOP (6) 

DOUBLE PRECISION :: APPLY (6), APPLY_E(6), C(3, 3), CN(3, 3), UN (3) 
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! — Determine 3x3 normal component stiffness matrix 


C(l, 

1) 

= E* ( 1 

- XNU ) / ( 1 - XNU - 2*XNU**2) 

C(l, 

2) 

= E*XNU/(1 - XNU - 2*XNU**2) 

C(l, 

3) 

= C(l, 

2) 

c (2, 

1) 

= C(l, 

2) 

C (2, 

2) 

= C(l, 

1) 

C (2, 

3) 

= C(l, 

2) 

C (3, 

1) 

= C(l, 

2) 

C (3, 

2) 

= C(l, 

2) 

C (3, 

3) 

= C(l, 

1) 


! — Rearrange equations such that all unknowns are on r.h.s. 
DO 230 I = 1, 3 

IF (LOP (I) .EQ. 2) THEN 

DO J = 1, 3 

IF (J .NE. I) CN ( I , J) = - C(I, J) / C(I, I) 
END DO 

CN ( I , I) = 1.0 / C (I, I) 

! — Cycle through other rows 
DO 210 K = 1, 3 

IF (K .NE. I) THEN 

DO J = 1, 3 

IF (J .EQ. I) THEN 


CN (K, 
ELSE 

J) 

= C (K, 

J) 

/ C(I, 

I) 

CN (K, 

J) 

= C (K, 

J) 

- C (K, 

I) * C(I, J) / C(I, I) 


ENDIF 
END DO 

ENDIF 

210 CONTINUE 

! — Reset C for next time (i.e. next specified S) 
DO II = 1, 3 
DO JJ = 1, 3 

ecu, jj) = cn(h, jj) 

END DO 
END DO 

ENDIF 
230 CONTINUE 


APPLY_E = APPLY 

IF (LOP ( 1 ) .EQ. 1) APPLY_E ( 1 ) = APPLY (1) - E11P 

IF (LOP (2) .EQ. 1) APPLY_E ( 2 ) = APPLY (2) - E22P 

IF (LOP (3) .EQ. 1) APPLY_E ( 3 ) = APPLY (3) - E33P 

! -- Solve for unknowns 
DO I = 1, 3 
UN (I) = 0 
DO J = 1, 3 

UN (I) = UN (I) + C(I, J) *APPLY_E ( J) 

END DO 
END DO 

IF (LOP ( 1 ) .EQ. 1) THEN 
EPS11 = APPLY (1) 

SIGH = UN ( 1 ) 

ELSE 

EPS11 = UN ( 1 ) + EllP 
SIGll = APPLY (1) 

ENDIF 

IF (LOP (2) .EQ. 1) THEN 
EPS22 = APPLY (2) 

SIG22 = UN (2) 

ELSE 

EPS22 = UN (2) + E22P 
SIG22 = APPLY (2) 

ENDIF 

IF (LOP (3) .EQ. 1) THEN 
EPS33 = APPLY (3) 

SIG33 = UN (3) 

ELSE 

EPS33 = UN (3) + E33P 
SIG33 = APPLY (3) 

ENDIF 

IF (LOP (4) .EQ. 1) THEN 
EPS23 = APPLY (4) 

SIG23 = 2*C44*APPLY_E (4) 

ELSE 

EPS23 = APPLY_E (4) / (2*C44) + E23P 
SIG23 = APPLY (4) 

ENDIF 
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IF (LOP (5) 
EPS13 = 
SIG13 = 
ELSE 

EPS13 = 
SIG13 = 
ENDIF 


.EQ. 1) THEN 
APPLY (5) 

2*C44*APPLY_E (5) 

APPLY_E (5) / (2*C44 ) 
APPLY (5) 


+ E13P 


IF (LOP (6) 
EPS12 = 
SIG12 = 
ELSE 

EPS12 = 
SIG12 = 
ENDIF 


.EQ. 1) THEN 
APPLY (6) 

2*C44*APPLY_E (6) 

APPLY_E (6) / (2*C44) 
APPLY (6) 


+ E12P 


! — Determine effective stress according to applied load 

SEFF_LOAD = DSQRT((SIG11 - SIG22)**2 + (SIG22 - SIG33) **2 + (SIGH - SIG33) **2 + & 
6 . * (SIG23**2 + SIG13**2 + SIG12**2) ) / DSQRT(2.D0) 


RETURN 

END 
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